Formation of defects in multirow Wigner crystals 
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We study the structural properties of a quasi-one-dimensional classical Wigner crystal, confined 
in the transverse direction by a parabolic potential. With increasing density, the one-dimensional 
crystal first splits into a zigzag crystal before progressively more rows appear. While up to four rows 
the ground state possesses a regular structure, five-row crystals exhibit defects in a certain density 
regime. We identify two phases with different types of defects. Furthermore, using a simplified 
model, we show that beyond nine rows no stable regular structures exist. 

PACS numbers: 61.50.-f, 71.10.Pm 



INTRODUCTION 



The electron crystal has created considerable inter- 
est since its possible existence was first pointed out by 
Wigner—. The three-dimensional Wigner crystal and its 
two-dimensional counterpart have been extensively stud- 
ied, and there exist beautiful experimental realizations of 
the latter using electrons trapped on the surface of liquid 
Helium^Ti^. More recently, Wigner crystallization in one 
dimension has received renewed interest^"—; for recent 
reviews see Refs. [2^ and [23 . 

The realization of a one-dimensional system requires 
the dominance of the confining potential over internal en- 
ergies, in particular, the inter-particle interactions. Upon 
increasing density (and, thus, the interaction energy), or 
weakening the confining potential, the crystal deviates 
from its strictly one-dimensional structure. It has been 
shown that at a critical density, a transition to a zigzag 
crystal takes placei'iii^'2^1^. Though not for electrons, 
this zigzag transition has indeed been observed using 
^^Mg"*" ions in a quadrupole storage ring22. 

Here we investigate the structural properties of the 
classical quasi-one-dimensional Wigner crystal beyond 
the zigzag regime. While previous investigations^ con- 
centrated on regular structures, we are interested in the 
formation of defects. From symmetry considerations the 
assumption of regular crystals is plausible at low densities 
when the number of rows is small, however, its validity is 
not at all obvious once the lateral extent of the crystal in- 
creases at higher densities. In fact, one expects a nonuni- 
form charge density in the direction transverse to the wire 
axis. In particular, considering the electrostatics problem 
of charges confined by a parabolic potential, V{y) oc y^, 
the density profile should obey n{y) cx -y/w^ — y^, where 
w is the width of the system^^'^^. Therefore, the assump- 
tion of perfect rows with equal linear densities should 
eventually break down. The formation of defects is of 
particular interest because they will have a direct impact 
on the transport properties of the system: while regular 
rows are locked, defects are expected to be mobile. 



II. MODEL 

We consider classical particles in two dimensions inter- 
acting via long-range Coulomb interaction. The system 
is assumed to be infinite in the x-direction and confined 
in the transverse y-direction by a parabolic confining po- 
tential Vconi- The energy of the system then reads 
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where e is the dielectric constant of the material and f2 
is the frequency of harmonic oscillations in the confining 
potential. 

At low densities, the system is one-dimensional, and 
the particles minimize their mutual Coulomb repulsion 
by occupying equidistant positions along the wire, form- 
ing a structure with short-range crystalline order-the so- 
called one-dimensional Wigner crystal"'^. Upon increas- 
ing the density, the inter-electron distance diminishes, 
and the resulting stronger electron repulsion eventually 
overcomes the confining potential V^onf , transforming the 
classical one-dimensional Wigner crystal into a staggered 
(zigzag) chain. From the comparison of the Coulomb in- 
teraction energy Vint('') = e^/er with the confining po- 
tential an important characteristic length scale emerges. 
Indeed, the transition from the one-dimensional Wigner 
crystal to the zigzag chain is expected to take place when 
distances between electrons are of the order of the scale 
To such that Vconf(''o) = Mnt(''o)- Within our model, i.e., 
for a parabolic confining potential and Coulomb interac- 
tions, the characteristic length scale tq is given as 
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It is convenient for the following discussion to measure 
lengths in units of tq. To that purpose we introduce a 
dimensionless density 

V = nero, (4) 
where rig = N/ L is the linear density of the system. 
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TABLE I: Number of rows in the crystal as a function of the 
dimensionless density v, assuming regular structures. 
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1.79 < !^ < 2.72 
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2.72 <v < 3.75 


5 


3.75 <u < 4.84 
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4.84 < !^ < 5.99 



Rescaling lengths, the energy can be written as 
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where Eq ^ {e^m^^ j^f?)"'^ . 

As a first step, we minimize the energy assuming regu- 
lar rows, aiming to find approximate values for the den- 
sity range in which a configuration with a given number 
of rows is stable. Assuming staggering in the x-direction 
between neighboring rows and inversion symmetry of the 
y-positions of the rows with respect to the wire axis, the 
number of minimization parameters is M/2 ((Af — l)/2) 
for even (odd) number of rows M, and the minimiza- 
tion is straightforward. Within these constraints, the 
minimization of the energy with respect to the electron 
configuration reveal o'^i^i^^ that a one-dimensional crystal 
is stable for densities v < 0.78, whereas a zigzag chain 
forms at intermediate densities 0.78 < ly < 1.71. More 
rows appear as the density further increases. The num- 
ber of rows as a function of v is shown in Table HI One 
notices that, with the exception of the four- row struc- 
turei in the regime 1.71 < u < 1.79, the linear density 
per row v^ow = v/M is of order < 1 in all cases, i.e., 
another row is added to the crystal when the distance 
between particles within a row is of the order of rg. A 
typical regular structure is shown in Fig. [TJ 

To investigate the importance of defects, the above 
conditions have to be relaxed. In the following, we con- 
centrate on the density regime 1.79 < < 5.99, encom- 
passing structures with 3 to 6 rows. In Sec. IIIII the nu- 
merical method is introduced, and in Sec. II VI we present 
our results. In Sec.|V]we introduce a simplified minimiza- 
tion procedure that allows us to extend the calculation 
to a larger number of rows, before concluding in Sec IVII 



III. NUMERICAL METHOD 

In order to find the ground state configuration of the 
system, the energy of the electrons in the parabolic con- 
fining potential is minimized with respect to the posi- 
tions of the electrons for given confinement strength and 




FIG. 1: Regular structure with 5 rows at v = 4.17, shown 
with its Voronoi construction for illustration purposes. This 
structure was obtained for 60 electrons in the unit cell. 



density. As the number of particles used in the simu- 
lation is finite, commensurability effects are important. 
To realize a regular Af-row structure, the number of par- 
ticles in the simulation box has to be a multiple of M. 
Similarly, to realize a defected structure, the defect den- 
sity is determined by the number of particles used in the 
simulation. To illustrate this, let us consider a five-row 
structure. Regular structures are realized for = 5n; 
for all other A^, defects appear. As we expect the density 
to be maximal at the center and decrease towards the 
edges, the simplest symmetric defected structure possi- 
ble is one where the outer rows are missing one particle 
each compared to the inner rows, i.e., structures of the 
form [{n—l)nnn{n — l)]. Such structures are realized 
for N — 5n — 2. The defect density may be defined 
as the number of missing particles in the outer rows 
divided by the number of particles in the inner rows, 

"dcf = ("inner " 'raouter)/ninnor = 1/n = 5/(A'' -|- 2). The 

minimum defect density that can be realized is, therefore, 
determined by the maximal number of particles that can 
be simulated. Thus, to find the ground state of the sys- 
tem, we have to vary N at fixed confinement strength 
and density. 

Conceptually, the proposed calculation is straightfor- 
ward. The computational difficulty arises from the com- 
plexity of the minimization problem. It is well known 
from the study of related problems, e.g., the determina- 
tion of the ground state of atomic clusters or the opti- 
mal arrangement of charges in a two-dimensional con- 
fined geometrjj^"— , that the corresponding energy func- 
tional has a number of metastable states that increases 
exponentially with the number of particles. In such a 
case, classic minimization techniques are not the optimal 
choice. 

Hybrid techniques employing genetic algorithms have 
been used in many related problems^ '"— as a general tool 
to explore the available phase space more thoroughly and 
obtain better solutions with comparable computational 
cost to conventional optimization techniques. One fre- 
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quently finds that counterintuitive disordered structures 
are favored. 

For our case, a simulation box of finite lengtlr L con- 
taining N electrons is used. Periodic boundary condi- 
tions in the x-direction are enforced to remove size ef- 
fects. For the summation of the interaction series, a 
quasi-one-dimensional restriction of the Ewald method is 
employed, following a similar technique to that reported 
in Ref. ISSl . The appropriate methods of proven stability 
for our quasi-one-dimensional geometry are of complexity 
O(iV^) and this fact, in conjunction with the significant 
number of minimizations that need to be carried out (var- 
ious system sizes for given total linear density), implies 
the necessity of substantial computational resources. 

The total energy per particle of a particular configura- 
tion of N electrons {r^ } can be written as 

where Eq is the previously defined energy scale and dis- 
tances are now measured in units of L. The complicated 
expression for the (dimensionless) interaction energy e 
and the details of its calculation are shown in appendix 
1X1 For a given number of electrons in a cell of length 
L, and a given density one has to minimize E [{ry }] 
with respect to the electron configuration and thereby 
obtain the stable structure with energy i?Gs(i^; 

In a nutshell, the algorithm proceeds along the follow- 
ing steps: An initial population of structures with ran- 
dom arrangements of electrons within the cell is partially 
relaxed towards a (local) minimum by a small number of 
iterations of a conventional minimization algorithm. Ev- 
ery member of the original population is then randomly 
split into two pieces, and the next generation is created 
by merging the pieces in all possible combinations while 
conserving the total number of particles. Subsequently, 
all newly obtained structures are fully relaxed to a (per- 
haps only local) minimum by a conventional minimiza- 
tion algorithm. A number of them is then chosen as 
parent structures for the next generation, always main- 
taining an appropriate diversity in the available configu- 
rations, i.e., a wide enough distribution in energies. The 
structure with the minimum energy is always retained 
to serve as a parent. The entire cycle is repeated un- 
til acceptable convergence is achieved. As expected, this 
hybrid approach is superior to simple minimization: it 
rapidly and consistently converges to complicated struc- 
tures, avoiding being trapped in local minima. 

In the end, to find the ground state configuration of 
the system at a given density the structure with the 
lowest energy, Eqs ~ minAr{£'Gs(i'; A'^)}, is chosen. 

IV. RESULTS 

With the method described above, we are able to con- 
sider systems comprised of up to ~ 200 electrons in the 



unit cell. We find that the lowest energy structures for a 
given energy are either regular structures, or structures 
where the linear density of the outer- most rows, I'outcr, is 
lower than the linear density of the inner rows, i^innei^- 
The finite number of particles in the unit cell implies a 
lower limit to the defect density we can consider. Here 
we define the defect density as njcf = 1 — t'outcr/i^innor- 
For up to 6 rows, the number of particles per row exceeds 
30. We are, therefore, able to identify defected structures 
with linear defect densities ridcf down to ~ 0.03. 

Let us summarize our main findings before discussing 
them in more detail: Up to 4 rows, the ground state of 
the system is free of defects. In the five-row structure, 
defects appear as one approaches the transition to 6 rows. 
Typical examples of such defected structures are shown 
in Fig. [21 We find that the defect density quickly in- 
creases with density and then levels off at values of the 
order ridcf ^ 0.08. Note that different types of defects ap- 
pear: In the low density regime where the defect density 
rapidly increases with density, the structure possesses in- 
version symmetry with respect to x-axis, i.e., the centers 
of the defects in the two outer rows are located at the 
same x-position as shown in Fig. [2^. By contrast, the 
structures with the maximal defect density ridof ~ 0.08 
display defects that are maximally shifted with respect 
to each other along the x-direction as shown in Fig. ^p. 
The transition to six-row structures is shifted to a larger 
density as compared to the value given in Table HI Above 
the transition, the ground state is a regular six-row struc- 
ture. Only upon further increasing the density do de- 
fects appear again, before the transition to a seven-row 
structure. Further analyzing the spatial structure of the 
ground state configurations, we find that the presence of 
defects in the outer rows also affects the particle posi- 
tions in the inner rows. While structures without defects 
consist of straight rows without corrugation, structures 
with defects display corrugation, i.e., distortions of the 
regular structure in both x- and j/-direction. 

A. Defects in Ave- and six-row crystals 

Using the full numerical minimization procedure, we 
find that the five-row Wigner crystal is stable in the den- 
sity range 3.75 < v < 4.86. Defected structures replace 
the regular ground state aX Vc — 4.695 and persist until 
the transition to 6 rows. Note that, as the finite number 
of particles limits the defect densities we can probe, the 

(5) 

value Vc represents an upper boundary for the range of 
stability of the regular structure. 

A regular five-row crystal is shown in Fig. [H whereas 
two defected five-row crystals are shown in Fig. |21 The 
five-row crystals in Fig. [21 correspond to = 4.70 and 
V = 4.85, and have a defect density ndcf — 0.038 and 
i-def = 0.083, respectively. Fig.Oshows how the defected 
structures were identified. In particular, the ground state 
energy at fixed density v = 4.85 as a function of the 
number of particles in the unit cell is shown. Regular 
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v=4.70 




a) 



FIG. 2: Defected structures with 5 rows at a) v = 4.70 and 
b) 4.85. The unit ceU consists of 128 and 58 electrons, re- 
spectively, with the two outer rows missing an electron. The 
corresponding defect densities are n^'™ = 1/26 ~ 0.038 and 
n^^f^ = 1/12 « 0.083. Electrons in red half-filled disks indi- 
cate the formal centers of the defects encountered. 



five-row structures are realized for N — 5n with rt G N. 
For all other N, defected structures are obtained. Let 
us discuss the high density structure v = 4.85 displayed 
in Fig. first. Three equivalent minima at = 58, 
N = 116, and N = 174 can be clearly seen. These 
minima correspond to configurations with defects in the 
outer- most rows, i.e., the outer rows have less particles 
than the inner rows, namely A^inncr — 12 (24,36), and 
A^outcr — 11 (22,33), respectively. The corresponding 
defect density is n^;f/ = 1 - 11/12 = 0.083. The low- 
energy structure v = 4.70 displayed in Fig. displays 
only one minimum at = 128 within the regime that 
can be explored by the full minimization procedure. This 
minimum corresponds to a defect density n'^'Jf^ = 0.0385. 
To rule out finite size effects, we extended our calcu- 
lation to a larger number of particles employing con- 
ventional minimization techniques, utilizing as starting 
guesses the structures obtained form the full minimiza- 
tion in the smaller unit cell. As Fig. [3^ shows, further 
minima appear at iV = 251, N = 379, and N = 502 cor- 
responding to approximately the same defect density. In 
fact, the lowest energy structure is obtained for N = 379 



where 



,4.70 _ 
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For comparison, Fig.|3]shows the equivalent diagram at 
a lower fixed density j/ = 4.17 corresponding to the regu- 
lar ground state shown in Fig. [TJ The energy of defected 
structures keeps decreasing with defect density until the 
lowest defect density ndcf = 1 — 39/40 — 0.025 reached 
given our limitation on the number of particles. Note 
that the lowest excitation branch shown in the picture 
corresponds to structures missing one particle from only 
one of the outer rows. Fitting that branch to a general 
functional form a + /3/N^ we obtain 7^1 and an energy 
gap in the thermodynamic limit AEoo — 2x 10~^. Thus, 



X 10 



v=4.70 




100 150 200 250 300 350 400 450 500 

N 



b) 



X 10 



v=4.85 



-1.5 



-2 - 



o o 



O o ' 



o o 
° • o ° 



50 



100 



150 



200 



N 



FIG. 3: Egs{N) at fixed density for a.) u = 4.70 and h) v = 
4.85. The energy is measured in units of iJo = {e'^mQ,'^ /2e'^)^^^ 
from the lowest-energy regular structure, a) For v = 4.70, 
four minima can be clearly seen. As the defect densities sam- 
pled differ slightly, they are not exactly equal in energy. The 
global minimum is found at A = 379, corresponding to a de- 
fect density n^e™ = 0.390. b) For v = 4.85, three degenerate 
minima appear at A = 58, 116, 174, all corresponding to a 
defect density of nifi" = 0.083. 



we do not expect that the regular structure becomes un- 
stable at even lower defect densities. 

Our findings in the vicinity of the transition from 5 to 6 
rows are summarized in Fig.[5l The defect density as well 
as the energy gaps to the lowest-lying regular or defected 
structure are shown as a function of density. Note that 
due to the substantial computational effort involved, the 
density interval is not uniformly sampled. As mentioned 
earlier, the defect density quickly increases in a narrow 
density interval and then levels off to an almost constant 
value until the transition to six rows is reached. The six- 
row crystal is stable in the density range 4.86 < i' < 6.04. 
It develops defects at around v — 5.75, which also persist 
until the transition to 7 rows. 

To better understand the structures that appear we 
now turn to a more detailed analysis of the spatial ar- 
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FIG. 4: Difference in energy at fixed density v — 4.17 as a 
function of the size of the unit cell. A'^ = 5n — 1 structures 
[{n — l)nnnn] form the lowest excitation branch, followed by 
iV = 5n + 1 [n{n + l)nnn], N = 5n - 2 [{n~l)nnn{n-l)], 
and Ai" = 5n — 3 [{n — l){n — l)nn(n — l)]. Fitting the lowest 
excitation branch to a general functional form a + jS/N"' we 
obtain 7 ~ 1 and an energy gap in the thermodynamic limit 
AE^ = 2 X 10"'^. 



rangement of particles in the crystal. 



B. Analysis of row corrugation 

As can be seen from Fig.[2l two types of defected struc- 
tures appear. These two structures can be distinguished 
by analyzing the distortion of the crystal. The distances 
between rows vary as a function of density. While reg- 
ular structures consist of straight rows, structures with 
defects display corrugation. Let us label the positions 
of particles as r^'^-', where j denotes the row and k de- 
notes the position along the row. In regular structures, 

we find yj*'' = yj*' for all j,k,k', within the accuracy 
of the calculation. For the defected structures, we define 
the average displacement of each row 



TV,- 



(7) 



fc=i 



where Nj is the number of particles in row j. In Fig. ^ 
the average positions of the rows are shown. Due to 
the symmetry of the structure only half of them are dis- 
played. 

As expected, there are jumps at the transition to a 
structure with a larger number of rows; in the interme- 
diate region the distance grows linearly with density. In- 
terestingly, as shown in the inset, a continuous transition 
to the structures with defects appears to take place. The 
transition also marks the onset of corrugation in the crys- 
tal structure. However, within the density regime of de- 
fected structures, we find a discontinuity. In the region 
labeled I, the distance between rows increases rapidly. 
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FIG. 5: Energy gaps |AiJ| and boundaries for the various 
phases encountered. The red vertical lines show the phase 
boundaries as obtained by the location of the zero of the cor- 
responding energy difference between: five-row regular and 
five-row defected structures (blue circles) , six-row regular and 
five-row defected structures (red rhombi) . The dashed line in- 
dicates the location of a first order transition between the two 
types of defects encountered: in region I the centers of the de- 
fects coincide whereas in region II the centers of the defects 
are maximally separated within the unit cell. The grey line 
shows the energy difference between five- and six-row regular 
structures. (Its zero is the prediction for the phase bound- 
ary under the assumption of regular rows.) Furthermore, the 
black stars show the defect density Tidof as a function of den- 
sity u. Note the jump at the boundary between regions I and 
II. 



a behavior that is well fitted by a square root. At the 
boundary between regions I and II, the row position dis- 
plays a jump before it increases linearly again in region 
II. This suggests two different defected phases that can 
be characterized by their corrugation. 

A close look at the defected structure reveals that the 
corrugation exists in both directions, along and perpen- 
dicular to the wire axis. We define the corrugation in the 
y-direction as the deviation from the average row position 
in that direction. 



(fc) _ 



(8) 



We can also define the average inter-particle distance for 
a given row j by 



(9) 



k=l 



where Vj is the dimensionless density in that row. Note 
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FIG. 6: Average positions of the crystal rows in units of ro as 
a function of density. Note that due to symmetry, only half of 
the rows are shown. The inset shows the detailed behavior at 
the transition from the regular five-row to the defected five- 
row structure. For regular structures, the distance between 
rows increases linearly with v. In region I, the distance be- 
tween rows increases more rapidly. At the boundary between 
regions I and II, there is a jump. In region II, the distance 
increases again linearly with the same slope as for the regular 
structures. 
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FIG. 7: Corrugation in the transverse (right panels) and longi- 
tudinal (left panels) direction for the defected structures that 
appear at a) v = 4.70 (top panels) and h) v = 4.85 (bottom 
panels). For the definition of the corresponding quantities 
please see the main text. The corrugation in the longitudi- 
nal direction (left) is measured in units of the average inter- 
particle distance within the row, while that in the transverse 
direction (right) is measured in units of the average distance 
between rows. The arrows indicate the particle located at the 
center of the defect in each case. See also Fig. [21 



that = X]n=i '^3- Subsequently, we define the corruga- 
tion in the cc-direction by 

6{Ax)f+' = xf+'^ - xf ^ - A^,. (10) 

While the corrugation is less than one percent in both 
directions, it turns out to be very important in deter- 
mining the ground state of the system. Figure [7] shows 
examples of the two dominant types of corrugation ac- 
companying five-row defected structures. Note that the 
arrows indicate the particle located at the center of the 
defect in each case. As before, the chosen density val- 
ues are v = 4.70 and v = 4.85, close to the boundaries 
shown in Figs. [S] and HI Fig. [7^ shows the corrugation for 
!/ = 4.70, i.e., a structure close to the density where de- 
fects first appear. This kind of corrugation is typical for 
the narrow density regime 4.695 < v < 4.712, where the 
defect density rapidly increases with i'. Fig.[7)D shows the 
corrugation for ly — 4.85 with defect density ndof — 0.083. 
This kind of corrugation is characteristic of the structures 
exhibiting the maximum defect density ridef 0.08. 

Qualitatively, the two types of structures exhibit dif- 
ferent features. In the first defected structure that is 
encountered, see Fig. [7^, the defects in the exterior rows 



are rather localized, and they are located at the same po- 
sition along the crystal. The displacements are maximal 
for the outer rows and decrease as one moves towards the 
interior of the crystal. In particular, due to the symme- 
try of the defect, the innermost row exhibits no corru- 
gation in the y-direction at all. At higher density, both 
the X- and y-corrugations are approximately sinusoidal. 
Furthermore, the defects on the two exterior rows are 
maximally separated, i.e., they are shifted by half a pe- 
riod, see Fig. [7}3. For this kind of structures, the center 
line possesses the maximum amplitude of y-oscillations. 
A possible explanation is that, while the interaction en- 
ergy (which drives the corrugation) is only sensitive to 
the relative corrugation, a deformation of the inner rows 
entails a smaller change in confining potential energy. 

We, thus, encounter two distinct phases with defects. 
Fig. [8^ shows the energy as a function of defect density 
for different densities close to the boundary between the 
two phases. Two minima corresponding to the different 
types of defected structures can be clearly identified. The 
position of one of the minima changes rapidly with den- 
sity. This minimum corresponds to the type of defect 
encountered in the low-density regime. The position of 
the other minimum barely shifts with density. This min- 
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imum corresponds to the sinusoidal defects encountered 
in the high-density regime. At low density, it describes 
a metastable state. However, its energy with respect 
to the other minimum decreases with density until, at 

(5) 

= 4.712, it eventually becomes the global minimum 
and, therefore, the ground state. Both the defect density 
(Fig. and the distance between rows (Fig. ^ display 
a discontinuity at the transition. The transition between 
the two defected phases is, thus, of first order. 

The nature of the transition from regular to defected 
structures is more difficult to identify as it requires going 
to very low defect densities. The fact, that with decreas- 
ing density, the defect density becomes lower and lower, 
until we reach the minimal value we can simulate, sug- 
gests, however, that the transition might be second order. 
In order to approach this transition, simplified models 
that allow one to simulate a larger number of particles 
are required. This models, then, may also be used in 
order to extend the calculation to larger number of rows. 



V. SIMPLIFIED MINIMIZATION 
PROCEDURES 

The full minimization procedure is computationally in- 
tensive which sets practical limits on the size of the unit 
cell one can simulate. That in turn imposes constraints 
on the defect density. Therefore, simplified models that 
allow us to simulate a larger number of particles are 
worth investigating to gain a better understanding. We 
start with the simplest model possible, compare with the 
results of the full simulation described above, and then 
discuss possible improvements. 

Up to six rows, we find that defected structures have 
less particles in the outer rows. In the previous section, 
we pointed out that these defected structures display cor- 
rugation. As a first approximation, one may neglect this 
corrugation and assume that all rows are straight and 
regular, i.e., Sy'^'^^ = (5(Ax)^'^'''^ '^^ — 0. Defects are in- 
corporated by allowing the linear densities in the inner 
and outer rows to differ-in particular, the two outer rows 
have less particles, noutcr < 'T-inncr- The density of defects 
is controlled by the parameter A — rtinncr/^-outcr, i-c, the 
density of defects is then given as ridof = 1 — A~^. In 
that case, for a fixed defect density, one has a minimiza- 
tion involving 2M — 1 parameters, namely the y-position 
of the rows and their relative shifts in the x-direction. 
Assuming that defects are located in the outer rows, the 
calculation can be further improved by "unfreezing" the 
x-positions of particles in the outer rows. This is the 
method we will use in the following. Given the much 
reduced parameter space, a conventional minimization 
procedure is sufficient here. 

The findings for five-row structures are shown in 
Fig. |8^. Note that the model captures correctly the ap- 
pearance of defects in the five-row structure, and it also 
predicts a regular ground state for the four-row crystal. 
Furthermore, the maximal defect density is reproduced: 



5 rows (full calculation) 




0.06 0.07 
"def 



b) 



5 rows (simplified calculation) 




FIG. 8: a) Difference AE — -Edcfcct — B'rGguiar for five-row 
structures as a function of defect density ridef within a simpli- 
fied calculation scheme that incorporates defects by allowing 
the hnear densities in the inner and outer rows to be different. 
Note that there is a region before the critical v — 4.86 of the 
5 — >■ 6 transition where the defected structures become stable, 
b) Energy of five-row defected structures at various densities 
as a function of defect density ridef using the full minimiza- 
tion procedure. The minima are shifted to zero for clarity of 
presentation. Note that the minimum present at small defect 
densities is not captured by the simplified method. 



For the five-row structure, we can see that the defect 
density ridct ^0.08 leading to the minimal energy is in 
agreement with the full minimization. 

Analyzing the results in more detail, (expected) dis- 
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crepancies are found. Due to the constraints imposed, 
the energies of defected structures are too high. For 
the five-row structures for example, the simphfied model 
finds that defects appear around v = 4.78 whereas the 
full minimization reveals that defected structures become 
the ground state configuration already at a lower den- 
sity V = 4.695. Furthermore, the simplified model does 
not capture the rise in defect density up to the maxi- 
mal value. In particular, the simplified model completely 
misses the low-density regime with symmetric defects. 
Comparing Figs. and b, the additional minimum at 
low-defect densities present in Fig. [5^ is clearly absent 
in Fig. ISIj. Thus, the simplified method only captures 
the high-density phase with sinusoidal defects. The rea- 
son is most likely that it corresponds to a fairly smooth 
corrugation and, therefore, is still present when one im- 
poses straight rows. By contrast, the minimum at low 
defect densities is associated with fairly sharp features in 
the corrugation profile and may, therefore, be suppressed 
by imposing straight regular rows. In particular, it is 
straightforward to verify that, for constant linear den- 
sity in the inner rows, the defects on the outer rows will 
be maximally separated. In order to obtain a defected 
structure where the centers of the defects coincide, a lon- 
gitudinal distortion of the inner rows is indispensable. 

To summarize, the simplified model correctly repro- 
duces the typical defect density-though it overestimates 
the energy of the defected structures which therefore are 
stable only in a reduced density interval. However, the 
model does not reproduce the low-density defected phase 
and, therefore, can not be used to explore the nature of 
the transition from regular to defected structures. The 
method may be used to study the stability of regular 
structures for structures with more rows. We find that 
upon further increasing the number of rows, the regime of 
densities where the ground state contains defects widens. 
Under the assumption that only the outer rows contain 
defects, regular structures disappear completely once the 
number of rows exceeds nine, as is evident from table |lll 
As the simplified method overestimates the energy of de- 
fected structures and misses the phase with symmetric 
defects, it is likely that regular crystals cease to be the 
ground state already for a smaller number of rows. 

As is also shown in Table |lTl the typical defect den- 
sity increases with the number of rows and also slightly 
varies with density for a given number of rows. Note that 
considering structures of the type [(n— . .n(n— 1)], 
the defect densities obtained can only take the discrete 
values ridef = 

Eventually, one expects that more complicated struc- 
tures will appear. The simple configuration we stud- 
ied is in competition with structures where defects ap- 
pear away from the edges, such as structures of the type 
[{n—l){n—l)n . . . n{n—l){n—l)], for example. Detailed cal- 
culations within these simplified models reveal that such 
structures are indeed competitors for the ground state, 
but up to 13 rows such a minimum is not realized. 

To approach the transition between regular and de- 



TABLE II: Number of rows and defect density in the crystal as 
a function of the dimensionless density All numbers shown 
were obtained using the simplified minimization procedure 
described in Sec. |V] 



#of 


total 


density range 


Wdcf 


rows 


density range 


with defects 




4 


2.72 <u < 3.75 


N/A 


N/A 


5 


3.75 <v < 4.86 


4.77 <v < 4.86 


0.083 


6 


4.86 <iy< 6.04 


5.75 <u < 6.04 


0.091-0.100 


7 


6.04 <u < 7.31 


6.76 <u < 7.31 


0.100-0.111 


8 


7.31 <v < 8.64 


7.80 <v < 8.64 


0.100-0.125 


9 


8.64 <u < 10.01 


8.87 <v < 10.01 


0.111-0.125 


10 


10.01 <iy < 11.37 


10.01 <iy < 11.37 


0.125-0.143 


11 


11.37 <iy < 12.77 


11.37 <iy < 12.77 


0.125-0.143 


12 


12.77 <u < 14.19 


12.77 <iy < 14.19 


0.143-0.167 


13 


14.19 <iy < 15.67 


14.19 <iy < 15.67 


0.143-0.167 



fected rows, we use a different trick. An unbiased search 
for the global minimum is numerically costly because a 
simple minimization may get stuck in a metastable min- 
imum. However, if the initial guess of the electron con- 
figuration is sufficiently close to the global minimum, a 
simple minimization will converge. Having identified the 
structure of defects in region I, one may feed such struc- 
tures into a simple minimization at lower densities. The 
results of such a procedure have been included in Figs. [S] 
and [6l There, structures with defect densities down to 
ndcf ~ 0.03 were obtained with the full minimization, 
whereas structures close to the phase boundary with 
lower defect structures were obtained with the method 
described here. The results suggest that the defect den- 
sity indeed vanishes at the transition which points to a 
second order phase transition. 



VI. CONCLUSION 

We study quasi-one-dimensional systems of classical 
particles interacting via long-range Coulomb interactions 
and confined by a parabolic potential in the transverse 
direction. The ground state configurations are multi-row 
Wigner crystals where the number of rows is controlled 
by the density (or the strength of the confining potential) . 
We find that defects that accommodate the density vari- 
ation in the transverse direction appear once the number 
of rows exceeds four. 

Defected structures have less particles in the outer 
than in the inner rows. The full numerical minimization 
for five rows reveals that two distinct types of defected 
phases exist. Upon increasing density, the regular struc- 
ture at low-densities is replaced by a structure with sym- 
metric defects, i.e., where the center of the defect on the 
two outer rows is located at the same x-position. As the 
number of particles that can be simulated sets a lower 
limit on the defect density that can obtained, the full 
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minimization allows one only to provide an upper limit 
for the density at the transition from regular to defected 
structures. We extend our calculations to lower densi- 
ties by using structures with the type of defect described 
above as the input for a simple minimization. The results 

(5) 

indicate that the defect density vanishes at z^c — 4.695 
and that the transition is of second order. To obtain 
symmetric defects, the longitudinal distortion of the in- 
ner rows, namely an increased density at the center of 
the defect, is crucial. Any analytical description of the 
transition would have to take into account this distortion. 

Upon further increasing density, the defect density 
rapidly increases. At a critical density, t'jjjj — 4.712, 
structures with a different type of defect corresponding 
to a sinusoidal distortion of the rows with a phase shift of 
half a period between the two outer rows become become 
the ground state. This second regime is characterized by 
a defect density that barely varies with density and ex- 
tends up to the transition to six rows. The transition 
between the two defected phases with different symme- 
tries is first order. 

Simplified models neglecting the corrugation of the 
rows only capture this second defected phase. Thus, 
these models do not allow one to further investigate the 
nature of the phase transition from regular to defected 
structures. However, as this second phase occupies most 
of the density interval, they may be used to study the sta- 
bility of regular structures upon increasing the number of 
rows. We find that beyond nine rows, no stable regular 
structures exist. Taking into account that the simplified 
model overestimates the energy of defected structures, 
we expect that stable regular structures may disappear 
even earlier. 



where the Fourier transform of f{x) is defined as 



m = f 



(A2) 



Let us consider the function f{x) = e^^^^'-^^ *. By com- 
pleting the square and carrying out the Fourier transform 
integration, we obtain the fundamental equation 



+ 00 

E 



,-(p+nL) 



(A3) 



where the reciprocal lattice vectors are given by C? = 
27rr7i/L with m = 0, ±1, ... . The following definition of 
the incomplete F function is extensively used and there- 
fore given here for reference: 



r(//, ux"^) 



,2m 



+ 00 



(A4) 



The system we are considering contains a basic cell of 
length L with N electrons. The spatial extent in the 
y-direction is limited by the confining potential. In the 
x-direction, we impose periodic replications of the basic 
cell to avoid edge effects. The interaction energy per cell 
of the system can be written 



i^j n 



nLS^\ 



3 



^lE- 



1 



InLicl ' 



(A5) 

where qi is the charge of particle i , = — Vj , and the 
index n runs over replicas of the unit cell. The artificial 
separation of the terms is for our convenience. We then 
introduce the notation 



Acknowledgments 



1 



nix 



(A6) 



We would like to acknowledge stimulating discussions 
with K. A. Matveev, Yu. V. Nazarov, and A. Melikyan. 
Part of the calculations were performed at the Ohio Su- 
percomputer Center thanks to a grant of computing time. 
This work was supported by the U.S. Department of En- 
ergy, Office of Science, under Contract No. DE-FG02- 
07ER46424. 



Appendix A: Ewald summation method for a 
quasi-one-dimensional geometry 

The method we use essentially follows the steps out- 
lined in Ref. [s^ It is based on the Poisson summation 
formula relating summations over direct and reciprocal 
space, 



+ 00 +00 

J: /(nL) = l F 



L ^ 

m— — oo 



L 



(Al) 



for |r| ^Oand$o-E„=.ol/I"ix|. 

In what follows we will split the summations in direct 
and reciprocal space. To cancel the divergencies appear- 
ing in the above sums, we will assume a uniform neutral- 
izing background charge. 

Using Eq. (|A4I) , we obtain the following representation. 



(A7) 



1 /■+°° 



Here we will introduce an artificial separation constant 
a which will control the splitting of the summation be- 
tween direct and reciprocal space. We then have $(r) — 
$>(r) -I- $<(r), where 



1 /•+°° 

erfc(Q!|r + nLx\) 



E' 



|r -t- nLxl 



(A8) 
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and 



Thus, 



$<(r)-^y/" dtt-i/2e-(^+"^)'*e-«'*. (A9) 
To evaluate $^(r), we use Eq. (jA3p yielding 



$<(r) = 1 Ve'G. Tdft-ig 



(AlO) 



While for G ^ the integration yields incomplete Bessel 
functional, 

j\tt-^e~^e-y'' = Ko , (All) 

the G — term (denoted Iq in the following) is divergent 
and has to be treated separately. Using the substitu- 
tion z = a^/t and expanding the second exponential, 
one finds 



In = — lim 



-t-oo 



dz . 



g2 +2^ (-1)™ 



2m^— m 



i G-s-0 /i 

The divergent contribution comes from m = 0, namely 



+ 00 



lim / dzz — — 7 + ln4a^— lim In G 



The rest of the sum can be evaluated to 



dzz^^^ '—iay)^"'z-"' 



E 

m— 1 ' 

+ CXD 



ml 



(A12) 



2^iayr- = -7 - In(a^y^) - r(0, aV)- 



Jo = -y (27 + lim In G^ + ln(?//4) + T (O, a^y') | , 



and 



cl><(r) = ±^e'G.^o(^,aV)+/o.(Al3) 

G/O ^ 



Splitting up $0 in the same way, we find 



&>(^\ - \^ erfc(a|n|L) 
^° ' ^ \n\L 



(A14) 



and 



1/2 



G#0 ^ 

-— (7 - ln4a^ + lim InG^l - 



At this stage we put everything together, e[{rij}] = 
1 Si^^j lilj^i^ij) + ^Qj^Oi and combining various terms 
we obtain the result for the interaction energy per cell of 
the system. 



m— 1 



~rr n 1 ■c"^' crfcfalri,- + nLxl) 1 

i,j In 



G#0 



G2 



„2„,2 



(A16) 



2L 



b + In + r(o, a'y-,)] - ^ E 



E 



7 - ln4a'^ + lim InG 

G-!-0 



where the notation X^n' implies that for n = there is tron in the simulation box can be cast as follows 
no self-interaction term in the summation. For a charge 

neutral system, the last term vanishes. For a system ^ r /2~ 

of electrons, as the one under consideration, a uniform £[{rij}] = -^/[{f ij}] -1= — N'' 

positive neutralizing background will exactly cancel the * ^ 



7 — In 



L 



divergent term limc-i-o In G . (A17) 

We define a dimensionless separation constant through 
a = a/L and introduce dimensionless coordinates. Sub- 
sequently, the dimensionless interaction energy per elec- with 
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f[{r^,}] 



E 



+00 



2 cos {2TTqXij ) Kq 

q=l 



2 9 
r ~2 2 

5^2 ' ^ Vi] 



ln5V,-r(0,52yf,).(A18) 
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